{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Check if all the packages are installed or not\n",
    "cond = \"Gadfly\" in keys(Pkg.installed()) &&\n",
    "\"Colors\" in keys(Pkg.installed()) &&\n",
    "\"ODEInterface\" in keys(Pkg.installed());\n",
    "@assert cond \"Please check if the following package(s) are installed:\\\n",
    "    Gadfly\\\n",
    "    Colors\\\n",
    "    ODEInterface\"\n",
    "\n",
    "# Load all the required packages\n",
    "using Gadfly\n",
    "using Colors\n",
    "using ODEInterface\n",
    "@ODEInterface.import_huge\n",
    "loadODESolvers();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "######################## Function for saving plots #################################\n",
    "# Input:\n",
    "# fileName = Name of the file where the plot is to be stored\n",
    "#            (with or without extension)\n",
    "# f_e = Array containing function evaluations as columns for each solver\n",
    "# err = Array containing erros as columns for each solver\n",
    "# solverNames = Array containing the names of solvers used in respective order\n",
    "# plotSize = size of the plot to be created\n",
    "# \n",
    "# Values have been tuned for a graph similar to the one in \n",
    "# Solving Ordinary Differential Equations I by\n",
    "# Hairer, Ernst, Nørsett, Syvert P., Wanner, Gerhard\n",
    "# page: 252\n",
    "###################################################################################\n",
    "function savePlotPNG(fileName,f_e,err,solverNames,\n",
    "    plotSize=[30cm,30cm])\n",
    "    \n",
    "    numOfLayers = length(solverNames);\n",
    "    \n",
    "    if !contains(fileName,\".\")\n",
    "        fileName = string(fileName,\".png\");\n",
    "    end\n",
    "    \n",
    "    plotColorsHex = [\"#4D4D4D\",\"#5DA5DA\",\"#FAA43A\",\"#60BD68\",\n",
    "        \"#F17CB0\",\"#B2912F\",\"#B276B2\", \"#DECF3F\",\"#F15854\"];\n",
    "    plotColors = [parse(Colorant,c) for c in plotColorsHex];\n",
    "    \n",
    "    majorFontSize = 24pt;\n",
    "    minorFontSize = 20pt;\n",
    "    pointSize = 5pt;\n",
    "    \n",
    "    myplot = plot(Scale.x_log10,Scale.y_log10,\n",
    "        Coord.cartesian(xflip=true),\n",
    "        Guide.manual_color_key(\"Legend\",solverNames,plotColorsHex[1:numOfLayers]),\n",
    "        Guide.xlabel(\"error\"),Guide.ylabel(\"#Function Evaluations\"),\n",
    "        Guide.xticks(ticks=[0:-3:-13;]),Guide.yticks(ticks=[3:1:5.1;]),\n",
    "        Theme(major_label_font_size=majorFontSize,panel_stroke=colorant\"black\",\n",
    "        minor_label_font_size=minorFontSize,key_title_font_size=majorFontSize,\n",
    "        key_label_font_size=minorFontSize,key_position=:top,key_max_columns=1));\n",
    "    \n",
    "    for i = 1:numOfLayers\n",
    "        push!(myplot,layer(x=err[:,i],y=f_e[:,i],Geom.point,Geom.path,\n",
    "        Theme(default_color=plotColors[i],default_point_size=pointSize)));\n",
    "    end\n",
    "    \n",
    "    draw(PNG(fileName,plotSize[1],plotSize[2]),myplot)\n",
    "    return nothing\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Number of subdivisions of the rope\n",
    "global n = 40;\n",
    "\n",
    "# Define the system of ODEs\n",
    "function rope(t,x,dx)\n",
    "    n2 = n*n; # n^2\n",
    "    n3by4 = convert(Int64,3*n/4); # 3*n/4\n",
    "    \n",
    "    # Force in x-direction\n",
    "    Fx = 0.4;\n",
    "    # Force in y-direction\n",
    "    Fy = cosh(4*t-2.5)^(-4);\n",
    "    \n",
    "    # Compute required matrices\n",
    "    c = -cos(x[1:n-1]-x[2:n]);\n",
    "    cDiag = [1;2*ones(n-2);3];\n",
    "    C = spdiagm((c,cDiag,c),(-1,0,1));\n",
    "    \n",
    "    d = -sin(x[1:n-1]-x[2:n]);\n",
    "    D = spdiagm((-d,d),(-1,1));\n",
    "    \n",
    "    # Compute the inhomogeneous term\n",
    "    v = -(n2+n/2-n*[1:n;]).*sin(x[1:n])-n2*sin(x[1:n])*Fx;\n",
    "    v[1:n3by4] = v[1:n3by4] + n2*cos(x[1:n3by4])*Fy; \n",
    "    \n",
    "    w = D*v+x[n+1:2*n].^2;\n",
    "    u = C\\w;\n",
    "    \n",
    "    # Write down the system\n",
    "    dx[1:n] = x[n+1:2*n];\n",
    "    dx[n+1:2*n] = C*v + D*u;\n",
    "    \n",
    "    return nothing\n",
    "end\n",
    "\n",
    "# Initial Conditions\n",
    "t0 = 0.0; T = 3.723; x0=zeros(2*n);\n",
    "\n",
    "# Compute the \"reference solution\"\n",
    "opt = OptionsODE(OPT_EPS => 1.11e-16,OPT_RHS_CALLMODE => RHS_CALL_INSITU,\n",
    "OPT_RTOL => 1e-16,OPT_ATOL=>1e-16);\n",
    "(t,x_ref,retcode,stats) = dop853(rope,t0, T, x0, opt);\n",
    "\n",
    "if retcode != 1\n",
    "    println(\"Reference solution failed\")\n",
    "else\n",
    "    # Initialization for the loop\n",
    "    # f_e = function evaluations\n",
    "    f_e = zeros(Int32,89,3);\n",
    "    # err = error for last step using infinity norm\n",
    "    err = zeros(Float64,89,3);\n",
    "\n",
    "    # solverNames = names of the solvers used for the plot\n",
    "    solverNames = [\"DOPRI5\",\"DOP853\",\"ODEX\"];\n",
    "\n",
    "    # Compute all the solutions\n",
    "    for i=0:88\n",
    "        \n",
    "        # Set up the tolerance\n",
    "        Tol = 10^(-3-i/8);\n",
    "        \n",
    "        # Set up solver options\n",
    "        opt = OptionsODE(OPT_EPS => 1.11e-16,OPT_RHS_CALLMODE => RHS_CALL_INSITU,\n",
    "        OPT_RTOL => Tol,OPT_ATOL => Tol);\n",
    "\n",
    "        # Solve using DOPRI5\n",
    "        (t,x,retcode,stats) = dopri5(rope,t0, T, x0, opt);\n",
    "        # Check if solver was successful\n",
    "        if retcode != 1\n",
    "            printFlag = false;\n",
    "            break;\n",
    "        end\n",
    "        f_e[i+1,1] = stats.vals[13];\n",
    "        err[i+1,1] = norm(x_accurate[1:n] - x[1:n],Inf);\n",
    "\n",
    "        # Solve using DOP853\n",
    "        (t,x,retcode,stats) = dop853(rope,t0, T, x0, opt);\n",
    "        # Check if solver was successful\n",
    "        if retcode != 1\n",
    "            printFlag = false;\n",
    "            break;\n",
    "        end\n",
    "        f_e[i+1,2] = stats.vals[13];\n",
    "        err[i+1,2] = norm(x_accurate[1:n] - x[1:n],Inf);\n",
    "\n",
    "        # Solve using ODEX\n",
    "        (t,x,retcode,stats) = odex(rope,t0, T, x0, opt);\n",
    "        # Check if solver was successful\n",
    "        if retcode != 1\n",
    "            printFlag = false;\n",
    "            break;\n",
    "        end\n",
    "        f_e[i+1,3] = stats.vals[13];\n",
    "        err[i+1,3] = norm(x_accurate[1:n] - x[1:n],Inf);\n",
    "    end\n",
    "\n",
    "    # Save the plot in PNG format\n",
    "    if printFlag\n",
    "        savePlotPNG(\"RopeConvTest\",f_e,err,solverNames);\n",
    "    else\n",
    "        println(\"Cannot generate plot due to solver failure\")\n",
    "    end\n",
    "end"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.2",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
